Solution-TP-3

1. Correction TP3 - Structures de Contrôle¶

Contenu

• 1. Correction TP3 - Structures de Contrôle
  • 1.1. Exo 1 : Structure de sélection / Expression ternaire
    • 1.1.1. Question 1
    • 1.1.2. Question 2
    • 1.1.3. Question 3
    • 1.1.4. Question 4
    • 1.1.5. Question 5
    • 1.1.6. Question 6
  • 1.2. Exo2 : Nombre de Froude
    • 1.2.1. Question 1
    • 1.2.2. Question 2
    • 1.2.3. Question 3
  • 1.3. Exo 3 : match/case - rayon hydraulique
    • 1.3.1. Question 1
    • 1.3.2. Question 2
    • 1.3.3. Question 3
  • 1.4. Exo 4 : for - data
    • 1.4.1. Question 1
    • 1.4.2. Question 2
    • 1.4.3. Question 3
  • 1.5. Exo5 : Matrice delta
    • 1.5.1. Description du code
    • 1.5.2. Réécriture de la boucle imbriquée j
  • 1.6. Exo 6 : Matrices [0],[1], [I]
  • 1.7. Exo 7 : $\pi$ (for/while)
    • 1.7.1. Question 1
    • 1.7.2. Question 2
    • 1.7.3. Question 3
    • 1.7.4. Question 4
    • 1.7.5. Question 5
  • 1.8. Exo 8 : Racines imbriquées
  • 1.9. Exo 9 Hauteur d'eau - formule de Manning

1.1. Exo 1 : Structure de sélection / Expression ternaire¶

1. Expliquer la commande suivante:

   signe = 1 if x > 0 else -1 if x < 0 else 0
   

2. Tester la commande pour différentes valeurs de x. Ajouter une commande qui permet de définir x et une autre qui affiche signe sous ce format : signe(valeur de x) = valeur de signe

Exemple : signe(-1.2) = -1, signe(1.2) = 1, signe(0) = 0

3. Rapporter le nom de l'erreur soulevée par la commande si on définit x = "1.2" (de type `str').

4. Réécrir la commande en plusieurs commandes d'une structure de sélection équivalente (if/elif/else)

5. Ajouter un test qui permet de vérifier le type de x et de n'accepter que les types int et float. Pour tout autre type, le code doit afficher par exemple : x = '1.2' n'est pas de type numérique

6. Si on définit signe = x/abs(x) (ou signe = abs(x)/x), quelles précautions supplémentaires doit-on prendre. Traduire en test et modifier le code.

1.1.1. Question 1¶

La commande est une expression ternaire (ou expression conditionnelle). Elle permet de determiner le signe d'un nombre entier ou réel.

1.1.2. Question 2¶

Il faut faire des tests sur le code

In [1]:

x = 0.1
signe = 1 if x > 0 else -1 if x < 0 else 0
print(f"signe({x}) = {signe}" )

signe(0.1) = 1

1.1.3. Question 3¶

pour x = "1.2", le code et soulève une erreur de type : TypeError

1.1.4. Question 4¶

Structure de selection équivalente if/elif/else

In [2]:

x = 0
if x > 0:
    signe = 1
elif (x < 0):
    signe = -1
else:
    signe = 0
print(f"signe({x}) = {signe}" )

signe(0) = 0

1.1.5. Question 5¶

Ajout d'un test qui permet de vérifier le type de n et de n'accepter que intet float

In [3]:

x = 1.2
# on peut utiliser type() ou bien isinstance()
# isinstance est recommandé

#if type(x) in (int,float):
if isinstance(x, (int,float)):
    if x > 0:
        signe = 1
    elif (x < 0):
        signe = -1
    else:
        signe = 0
    print(f"signe({x}) = {signe}")
else:
    signe = None
    print(f"{x = } n'est pas de type numérique")

signe(1.2) = 1

In [4]:

# Même que précédent avec expression ternaire
x = '1.2'
if isinstance(x, (int,float)):
    signe = 1 if x > 0 else -1 if x < 0 else 0
    print(f"signe({x}) = {signe}")
else:
    print(f"{x = } n'est pas de type numérique")

x = '1.2' n'est pas de type numérique

1.1.6. Question 6¶

Si on utilise signe = x/abs(x) ou signe = abs(x)/x Il faut faire attention à x = 0

In [5]:

# Utilisation de abs(x)
x = 0
if isinstance(x, (int,float)):
    signe = 0 if x == 0 else x/abs(x)
    print(f"signe({x}) = {signe}")
else:
    print(f"{x = } n'est pas de type numérique")

signe(0) = 0

Remarque
On peut exploiter la fait que True / False sont interprétés comme 1 / 0, pour calculer signe avec l'expression : signe = (x>0) - (x<0).

• Si x>0 : alors (x>0) - (x<0) donnera 1 - 0
• Si x<0 : alors (x>0) - (x<0) donnera 0 - 1
• Si x==0: alors (x>0) - (x<0) donnera 0 - 0

In [6]:

x = 5
if isinstance(x, (int,float)):
    signe = (x>0)-(x<0)
    print(f"signe({x}) = {signe}")
else:
    print(f"{x = } n'est pas de type numérique")

# Ou bien, en plus compact :
signe = (x>0)-(x<0) if isinstance(x, (int,float)) else None
print(f"signe({x}) = {signe}")

signe(5) = 1
signe(5) = 1

1.2. Exo2 : Nombre de Froude¶

Le code sivant permet de classer le régime d'un écoulement à surface libre selon la valeur du nombre de Froude $Fr = V /\sqrt{g h}$ avec $g=9.81 \ m/s^2$.

L'écoulement est caractérisé par la vitesse $V$ et la hauteur d'eau $h$.

1. Tester différentes valeurs de $V$ et de $h$ de sorte à avoir les 3 cas du régime d'écoulement. Tester particuilièrement les 3 cas suivants puis commenter les résultats et tirer une conclusion

   • a) $V=2.18 \ m/s$, $h = 0.48 \ m$
   • b) $V=2.52 \ m/s$, $h = 0.65 \ m$.
   • c) $V=1.918 \ m/s$, $h = 0.375 \ m$.

2. Réécrire la structure de sélection (if-elif-else) en utilisant une expression ternaire (conditionnelle)

3. Améliorer le code en introduisant une tolérance dans les conditions sur Fr. Refaire les tests des cas précédents et faire d'autres tests en fonction du seuil de tolérance.

1.2.1. Question 1¶

Exemple de réponses :

cas 1

V, h = 1.8, 0.5
fr = 0.8
regime = fluvial

cas 2

V, h = 2.5, 0.7
fr = 1.26
regime = torrentiel

Conclusion :

Le cas de régime critique est presque impossible à obtenir ni à détecter à cause des erreurs numériques. Le teste d'égalité est trop stricte et dépond aussi des erreurs machine.

In [7]:

# Données : tester différentes valeurs de V et h
V = 2.5     # vitesse moyenne (m/s)
h = 0.4     # hauteur d'eau (m)
g = 9.81     # accélération gravitationnelle (m/s²)

Fr = V/(g * h)**(0.5) # calcul de Fr

if Fr < 1.0:
    regime = "fluvial"
elif Fr == 1.0:
    regime = "critique"
else:
    regime = "torrentiel"

print(f"{V = :.2f}, {h = :.2f}")
print(f"Régime : {regime}  ({Fr = :.2f})")

V = 2.50, h = 0.40
Régime : torrentiel  (Fr = 1.26)

1.2.2. Question 2¶

In [8]:

# Prise en compte d'une tolérance

V = 1.51     # vitesse moyenne (m/s)
h = 0.25      # hauteur d'eau (m)
g = 9.81     # accélération gravitationnelle (m/s²)

Fr = V/(g * h)**(0.5) # calcul de Fr

tol = 5e-2        # tolérance
if Fr < 1.0 - tol:
    regime = "fluvial"
elif Fr > 1.0 + tol:
    regime = "torrentiel"
else:
    regime = "critique"

print(f"{V = :.2f}, {h = :.2f}")
print(f"Régime : {regime}  ({Fr = :.2f})")

V = 1.51, h = 0.25
Régime : critique  (Fr = 0.96)

1.2.3. Question 3¶

Expression ternaire avec et sans tolérence

Avec un tolérance de tol=0.05, on detecte bien le cas du régime critique

In [9]:

#Expression ternaire
V = 1.91     # vitesse moyenne (m/s)
h = 0.375      # hauteur d'eau (m)
g = 9.81     # accélération gravitationnelle (m/s²)

#h = (V**2)/g
#V = (g*h)**0.5
Fr = V/(g * h)**(0.5) # calcul de Fr

tol = 5e-2
regime = "fluvial" if Fr < 1.0-tol else "torrentiel" if Fr > 1.0+tol else "critique"

# sans tolérence
#regime = "fluvial" if Fr < 1.0 else "critique" if Fr == 1.0 else "torrentiel"

print(f"{V = :.2f}, {h = :.2f}")
print(f"Régime : {regime}  ({Fr = :.2f})")

V = 1.91, h = 0.38
Régime : critique  (Fr = 1.00)

1.3. Exo 3 : match/case - rayon hydraulique¶

On reprend dans cet exercise l'exemple du calcul du rayon hydraulique de la section 1.2.4 du cours des structures de contrôle.

Travail à réaliser :

1. Reprendre le code et ajouter le cas d'une section trapézoïdale.

2. Remplacer la structure match/case par une structure if/elif/else

3. Proposer une autre organisation des données et écrire le code correspondant qui permet de réaliser le même tâche.

In [10]:

# Définition des sections dans un dict
sections = {
    'r' : {'type' : 'rectangulaire', 'b' : 2.0, 'h' : 1.0},
    't' : {'type' : 'triangulaire' , 'h' : 1.0, 'm' : 1.0},
    'c' : {'type' : 'circulaire'   , 'd' : 1.2},
    }

# Choix du type de section et tests de correspondance.
section = 'r'    # type de section ('r', 't' ou 'c')

# bloc de structure de sélection match/case
match section:
    case 'r':                      # cas de section rectangulaire
        b = sections['r']['b']     # largeur (m)
        h = sections['r']['h']     # hauteur d'eau (m)
        A = b * h
        P = b + 2 * h
        Rh = A / P
        print(f"Section {sections['r']['type']} : Rh = {Rh:.3f} m")

    case 'c':                    # cas de section circulaire
        pi = 3.14159265
        D = sections['c']['d']   # diamètre (m)
        A = pi * (D**2) / 4
        P = pi * D
        Rh = A / P
        print(f"Section {sections['c']['type']} : Rh = {Rh:.3f} m")

    case 't':                       # cas de section triangulaire
        h = sections['t']['h']      # hauteur d'eau
        m = sections['t']['m']      # talus
        A = m * h**2
        P = 2 * h * (1 + m**2)**(1/2)
        Rh = A / P
        print(f"Section {sections['t']['type']} : Rh = {Rh:.3f} m")

    case _:                        # cas par défaut (autres cas)
        print("Type de section inconnu.")

Section rectangulaire : Rh = 0.500 m

1.3.1. Question 1¶

Section trapézoidale

In [11]:

# Définition des sections dans un dict
sections = {
    'r' : {'type' : 'rectangulaire', 'b' : 2.0, 'h' : 1.0},
    't' : {'type' : 'triangulaire' , 'h' : 1.0, 'm' : 1.0},
    'c' : {'type' : 'circulaire'   , 'd' : 1.2},
    'z' : {'type' : 'trapézoïdale' , 'b' : 1.4, 'h' : 0.8, 'm' : 1.0}

    }

# Choix du type de section et tests de correspondance.
section = 'z'
match section:
    case 'r':                      # cas de section rectangulaire
        b = sections['r']['b']     # largeur (m)
        h = sections['r']['h']     # hauteur d'eau (m)
        A = b * h
        P = b + 2 * h
        Rh = A / P
        print(f"Section {sections['r']['type']} : Rh = {Rh:.3f} m")

    case 'c':                    # cas de section circulaire
        pi = 3.14159265
        D = sections['c']['d']   # diamètre (m)
        A = pi * (D**2) / 4
        P = pi * D
        Rh = A / P
        print(f"Section {sections['c']['type']} : Rh = {Rh:.3f} m")

    case 't':                       # cas de section triangulaire
        h = sections['t']['h']      # hauteur d'eau
        m = sections['t']['m']      # talus
        A = m * h**2
        P = 2 * h * (1 + m**2)**(1/2)
        Rh = A / P
        print(f"Section {sections['t']['type']} : Rh = {Rh:.3f} m")

    case 'z':                       # cas de section trapézoidale
        h = sections['z']['h']      # hauteur d'eau
        b = sections['z']['b']      # largeur à la base
        m = sections['z']['m']      # talus
        A = (b + m * h ) * h
        P = b + 2 * h * (1 + m**2)**(1/2)
        Rh = A / P
        print(f"Section {sections['z']['type']} : Rh = {Rh:.3f} m")

    case _:                        # cas par défaut (autres cas)
        print("Type de section inconnu.")

Section trapézoïdale : Rh = 0.481 m

1.3.2. Question 2¶

Structure if/elif/else

In [12]:

# Définition des sections dans un dict
sections = {
    'r' : {'type' : 'rectangulaire', 'b' : 2.0, 'h' : 1.0},
    't' : {'type' : 'triangulaire' , 'h' : 1.0, 'm' : 1.0},
    'c' : {'type' : 'circulaire'   , 'd' : 1.2},
    'z' : {'type' : 'trapézoïdale' , 'b' : 1.4, 'h' : 0.8, 'm' : 1.0}

    }

# Choix du type de section et tests de correspondance.
section = 'r'

if section == 'r':                 # cas de section rectangulaire
    b = sections['r']['b']     # largeur (m)
    h = sections['r']['h']     # hauteur d'eau (m)
    A = b * h
    P = b + 2 * h
    Rh = A / P
    print(f"Section {sections['r']['type']} : Rh = {Rh:.3f} m")

elif section == 'c':               # cas de section circulaire
    pi = 3.14159265
    D = sections['c']['d']         # diamètre (m)
    A = pi * (D**2) / 4
    P = pi * D
    Rh = A / P
    print(f"Section {sections['c']['type']} : Rh = {Rh:.3f} m")

elif section == 't':            # cas de section triangulaire
    h = sections['t']['h']      # hauteur d'eau
    m = sections['t']['m']      # talus
    A = m * h**2
    P = 2 * h * (1 + m**2)**(1/2)
    Rh = A / P
    print(f"Section {sections['t']['type']} : Rh = {Rh:.3f} m")

elif section == 'z':            # cas de section trapézoidale
    h = sections['z']['h']      # hauteur d'eau
    b = sections['z']['b']      # largeur à la base
    m = sections['z']['m']      # talus
    A = (b + m * h ) * h
    P = b + 2 * h * (1 + m**2)**(1/2)
    Rh = A / P
    print(f"Section {sections['z']['type']} : Rh = {Rh:.3f} m")

else:                           # cas par défaut (autres cas)
    print("Type de section inconnu.")

Section rectangulaire : Rh = 0.500 m

1.3.3. Question 3¶

Ouverte aux choix des étudiants. Tout est accepté pourvu que le code marche.

1.4. Exo 4 : for - data¶

On considère la série de données suivante :

data = [56.96, 98.96,  83.81, 73.55, 56.31, 81.88,
        75.93, 50.75, 103.04, 85.38, 54.28, 73.13,
        81.50, 90.17, 102.21, 82.77, 38.77, 75.37]

Ecrire un code qui permet de calculer la moyenne et l'écart-type de la série data en 3 versions suivantes:

Écrire un code Python permettant de calculer la moyenne et l'écart-type de la série data selon les trois versions suivantes :

1. Version 1 - parcours par valeurs :
   On utilise une boucle de type for val in data

2. Version 2 - prcours par indices :
   On utilise une boucle for de type for i in range(n): avec n la longueur de la série.

3. Version 3 - version compacte :
   On utilise la fonction sum() et une compréhension de liste (ou un générateur).

In [13]:

data = [56.96, 98.96, 83.81, 73.55, 56.31,  81.88, 75.93, 50.75, 103.04,
        85.38, 54.28, 73.13, 81.50, 90.17, 102.21, 82.77, 38.77,  75.37]

1.4.1. Question 1¶

In [14]:

# Version 1 : parcours par valeurs
n = len(data)
mu = 0
for i in range(n):
    mu += data[i]
mu /= n
sig = 0
for i in range(n):
    sig += (data[i] - mu)**2
sig = (sig/n)**0.5
print('Version 1 : parcours par valeurs')
print(f"{mu  = :.2f}")
print(f"{sig = :.2f}")

Version 1 : parcours par valeurs
mu  = 75.82
sig = 17.74

1.4.2. Question 2¶

In [15]:

# Version 2 : parcours par indices
mu = 0
for val in data:
    mu += val
mu /= n

sig = 0
for val in data:
    sig += (val - mu)**2
sig = (sig/n)**0.5

print('version 2 : parcours par indices')
print(f"{mu  = :.2f}")
print(f"{sig = :.2f}")

version 2 : parcours par indices
mu  = 75.82
sig = 17.74

1.4.3. Question 3¶

In [16]:

# Version 3 - compacte -
mu = sum(data)/n
sig = (sum((x - mu)**2 for x in data)/n)**0.5

print('Version 3 : compacte ')
print(f"{mu  = :.2f}")
print(f"{sig = :.2f}")

Version 3 : compacte 
mu  = 75.82
sig = 17.74

In [17]:

print( sum(data[i] for i in range(n))/n )
print( sum(val for val  in data)/n )

75.82055555555556
75.82055555555556

1.5. Exo5 : Matrice delta¶

1. Décrire ce que fait le code ci-dessous tout en précisant ce que représente la variable delta et comment elle est construite. ( Noter bien l'expression int(i==j) )

2. Réécrire la boucle sur j avec une compréhension de liste

In [18]:

n = 3
delta = []
for i in range(n):
    ligne = []
    for j in range(n):
        ligne.append(int(i==j))   # liste pour la ligne i
    delta.append(ligne)           # ajout de la sous-liste
print(f'{delta = }')

delta = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

1.5.1. Description du code¶

Le code permet de construire une liste de listes avec deux boucles imbriquées. Le nombre de listes est n et chaque sous liste est de même taille n. Ce qui représente une matrice.

Puisque chaque sous liste prend la valeur int(i==j) on obtient la matrice identité. delta représente donc la matrice identité (ou le symbole de Kronecker pour n=3)

Remarque :
Les éléments delta[i][j] peuvent être obtenus de plusieurs manières, en voici trois :

1. Transtypage (casting) explicite en int : delta[i][j] = [int(i==j)]

2. Expression conditionnelle : delta[i][j] = [1 if i == j else 0]

3. Expression booléenne avec casting implicite : delta[i][j] = [1 - ((i < j) + (i > j))]

   • Ici, les comparaisons i < j et i > j renvoient des valeurs booléennes (True ou False), l'opérateur + les convertit implicitement en entiers (1 ou 0) lors de l'addition.

   • Les parenthèses sont importantes à cause de la priorité des opérateurs :+ a une priorité plus haute que < et >.

   Sans les parenthèses (i<j + i>j) est le même que i < (j + i) > j qui correspond à une comparaison enchaînée : i < (j+i) and (j+i) > j. Cela renvoie toujours un booléen (True ou False) et non un entier.

   • On peut également réécrire l’expression sous forme : delta[i][j] = [1 - (i<j or i>j) ]

1.5.2. Réécriture de la boucle imbriquée j¶

Modification de la boucle sur j en expression de compréhension de liste

In [19]:

n = 3
delta = []
for i in range(n):
    delta.append([int(i==j) for j in range(n)])
print(f'{delta = }')

delta = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

Remarque :
Il est possible de remplacer la fonction append() par une concaténation de liste (opérateur +). Pour tout indice i, on ajoute la compréhension de liste [int(i==j) for j in range(n)] comme sous-liste, en utilisant une double paire de crochets [[...]].

In [ ]:

# Version avec delta +=[ [...] ]
n = 3
delta = [] 
for i in range(n):
    delta+=[[int(i==j) for j in range(n)]]
print(f'{delta = }')

delta = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
delta = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

Si l'on veut une double boucle for, il faut d'abord ajouter une sous-liste vide à chaque indice i : delta += [[]]. Ensuite, dans la boucle imbriquée d'indice j, on ajoute les éléments un à un à la sous-liste delta[i] en les mettant entre crochets pour former une liste contenant un seul élément.

In [66]:

# double boucle
delta = []
for i in range(n):
    delta +=[[]]    # ajoute une sous-liste vide
    for j in range(n):
        delta[i] += [int(i==j)]
        #delta[i] += [1 if i==j else 0]
        #delta[i] += [1 - ( (i<j) + (i>j) ) ]
        #delta[i] += [ 1 - (i<j or i>j)  ]
        
        
[print(str(l)) for l in delta]
#print(f'{delta = }')

print(True or False)

[1, 0, 0]
[0, 1, 0]
[0, 0, 1]
True

1.6. Exo 6 : Matrices [0],[1], [I] ↑¶

Le code suivant permet de construire 3 matrices en utilisant des listes de listes :

1. Récrire le code pour construire les 3 matrices en utilisant une seule boucle for
2. Améliorer l'affichage

In [20]:

m, n = 3, 3 # dimensions des matrices

# zeros : matrice de zeros
zeros =  [[0 for k in range(m)] for _ in range(n)]

# ones : matrice de uns
ones = [[1 for _ in range(m)] for _ in range(n)]

# identity : matrice identité
identity = [[1 if i == j else 0 for j in range(3)] for i in range(3)]


print('ones     =' , ones)
print('identity =' , identity)
print('zero     =' , zeros)

ones     = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]
identity = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
zero     = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]

In [21]:

# Réponse
I = []   # l'Identité
U = []   # les Uns
Z = []   # les Zeros
for i in range(n):
    ligneI = []
    ligneU = []
    ligneZ = []
    for j in range(m):
        if i==j:
            ligneI.append(1)
        else:
            ligneI.append(0)
        # ligneI.append(1) if i==j else ligneI.append(0) # Expression conditionnelle
        ligneU.append(1)
        ligneZ.append(0)
    I.append(ligneI)
    U.append(ligneU)
    Z.append(ligneZ)

print('les Uns    =' , U)
print("l'Identité =" , I)
print('les zeros  =' , Z)

les Uns    = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]
l'Identité = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
les zeros  = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]

1.7. Exo 7 : $\pi$ (for/while)¶

La commande suivante

pi = 4*sum( (-1 if n % 2 else 1) / (2*n + 1) for n in range(N))

permet de calculer approximativement la valeur de $\pi$ en utilisant $N$ termes de la série de Leibnitz : $$ \pi \approx 4\ \frac{(-1)^n}{2n+1} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} \cdots $$

1. Expliquer cette commande en précisant les points suivants :

• Le type de l'objet sur lequel est appliquée la somme sum()
• Le type de l'expression dans le numérateur et le résultat attendu
• L'expression du dénominateur et son résultat
• La boucle for et le range
• Pourquoi la commande a évité (-1)**n

2. Ajouter la commande qui définit N puis celle qui affiche pi avec 8 chiffres décimaux.

3. Executer le code obtenu pour différentes valeurs de N et rapporter (sous forme d'un tableau) les valeurs de pi obtenues en fonction de N

4. Réécrire le code en utilisant la boucle for tout en remplaçant le terme (-1 if n % 2 else 1) par autre équivalent différent de (-1)**n

5. Réécrire le code en utilisant la boucle while basée sur un critère de précision eps et non pas sur le nombre d'itérations N. Ajouter un conteur d'itérations it et noter sa valeur en fin de boucle pour différentes valeurs de eps (de 1e-5 à 1e-8)

1.7.1. Question 1¶

• La commande applique la fonction sum() (somme de python) sur une expression génératrice puisqu'on a une boucle for dans une expression entre parenthèses. Le générateur génère des nombres un à un à la demande sans mobiliser trop de mémoire systeme

• Le numérateur est une expression conditionnelle (ternaire) qui prend -1 si le reste de la vision de n par 2 n'est pas nul sinon 1. Ce qui donne -1 si n est impaire et +1 si n est paire. Cela reproduit les altérnances des signes (+) (-) dans la série

• Le dénominateur est une expression régulière simple qui produit les nombres impaires 1,3, 5, ....

• Le range(N) produit tous les entiers entre 0 et N

• La commande a évité (-1)**n car l'opérateur ** effectue une élévation à la puissance, qui est, elle seule, beaucoup plus coûteuse en temps de clcul que les deux commandes modulo (%) et test if/else réunies.

1.7.2. Question 2¶

Le code avec les 3 commandes :

In [22]:

N = 1_000_000
pi = 4*sum( (-1 if n % 2 else 1) / (2*n + 1) for n in range(N))
print(f'Valeur approximative de pi : {pi:0.8f}')

# En meilleur format
print(f'\u03C0 \u2248 {pi :.8f}')

Valeur approximative de pi : 3.14159165
π ≈ 3.14159165

Remarque :
On rappelle que $n$ modulo $2$ (n % 2) donne 0 si $n$ est pair et 1 si $n$ est impair. Donc (-1 if n % 2 else 1) donnera -1 si $n$ est impaire et 1 si $n$ est pair. Ce qui correspond bien au terme d'alternance $(-1)^n$.

On peut aussi reproduire l'alternance avec : 1-2*(n%2) :

In [23]:

# une autre façon de calculer (-1)^2
N = 1_000_000
pi = 4*sum( (1 - 2*(n%2) ) / (2*n + 1) for n in range(N))
print(f'Valeur approximative de pi : {pi:0.8f}')

Valeur approximative de pi : 3.14159165

1.7.3. Question 3¶

Quelques valeurs de pi en fonction de N :

---

N	valeur calculée	temps d'execution
1 000	3.14059265	0.001
50 000	3.14157265	0.009
100 000	3.14158265	0.015
500 000	3.14159065	0.075
1 000 000	3.14159165	0.143
10 000 000	3.14159255	1.463
100 000 000	3.14159264	15.451
200 000 000	3.14159265	31.562

Remarque 1 : Le temps d'execution n'est pas demandé, on le donne ici juste à titre indicatif. Il serait intéressant de comparer entre les commandes : 2%n et (-1)**n.

Remarque 2 : Le temps de calcul dépend de la machine et du système. Il varie aussi légèrement d'une execution à une autre selon les autres processus en cours d'execution

1.7.4. Question 4¶

Réécriture du code avec la boucle for

La somme est initialisé à zéro, et le numérateur ainsi que le dénominateur à 1.

A chaque pas de calcul, le rapport numérateur/dénominateur est caculé et cumulé dans la somme. Ensuite, le numérateur change de signe et le dénominateur est augmenté de 2.

Le fait d'utiliser le changement de signe du numérateur, on évite (-1)**n et on gagne en temps d'excution. Par exemple pour N = 100 000 000, la commande somme += (-1)**n /(2*n+1) induit un ralangement du temps d'execution de 4.7s à 13.s soit : 8.6 s qui correspond approximativement = 183%.

In [24]:

N = 2_000_000

numer = denom = 1  # initialisation de la fraction
somme = 0
for n in range(N+1):
    somme += numer/denom # cumul de la somme
    numer = -numer       # changement de signe
    denom += 2           # incrémentation avec un pas 2

print(f'{n = :_d}')      # Affichage
print(f'pi \u2248 {4*somme :.8f}')

from math import pi
print(f'{pi = :0.8f}')

n = 2_000_000
pi ≈ 3.14159315
pi = 3.14159265

1.7.5. Question 5¶

Réécriture du code avec la boucle while.

Pour contrôler la boucle while, il nous faut une condition qui utilise la précision imposée. Cette condition sera appaliquée sur les termes de la série. Tant que la valeur absolue du terme est supérieur à la précision eps la boucle reste active.

Comarativement à la boucke for, il nous faut calculer le 1er terme avant la boucle pour lui appliquer la condition abs(terme) > eps et enclencher les itérations. A l'intérieur, le reste ne change pas beaucoup.

Le compteur it permet de compter le nombre d'itérations effectuée. Il est facultatif mais très efficace pour arrêter la boucle en cas de non-convergence (cas où les termes de la série ne s'amenuisent pas).

On peut donc ajouter un test de dépassement d'un nombre important (1 milliard par exemple) d'itérations, auquel cas on casse la boucle avec break

La comparaison entre les temps d'execution des 3 versions des codes, sum-générateur et boucles for et while montre que la boucle for induit 31% de temps en plus et while induit 97% si on n'utilise pas break contre 174 lorsqu'on utilise break

In [25]:

# Calcul avec la boucle while.

eps = 5e-7               # Précision
numer = denom = 1
terme = numer/denom      # terme de la série

somme = 0                # initialisation de la somme
it = 0                   # Conteur d'iteration (facultatif)

while abs(terme) > eps:
    somme += terme       # cumul de la somme avec le terme
    numer = -numer       # changement de signe
    denom += 2           # incrémentation avec un pas 2
    terme = numer/denom  # mis à jour du terme de la série

    #it += 1                 # mis à jour de conteur
    #if (it >1000_000_000):  # facultatif mais sécuritaire
    #    break               # pour éviter la boucle infinie

print(f'{it = :_d}')
print(f'pi \u2248 {4*somme :.8f}')

#validation (c'est facultatif)
#from math import pi
#print(f'{pi = :0.8f}')

it = 0
pi ≈ 3.14159165

1.8. Exo 8 : Racines imbriquées¶

On reprend l'exercice 2.c du TP1 en le généralisant pour nombre réel $x$ et un nombre entier $k$ de termes sous racines

$r = \sqrt{x + \sqrt{x + \sqrt{x +\ ... }}} $

1. Ecrire un code qui utilise une boucle while pour calculer $r$ avec une precision eps donnée (prendre eps = 1e-16).

   Ce code doit permettre :

   • d'éviter le cas d'une boucle infinie avec un nombre d'iterations maximum $k_{max}$ à ne pas dépasser

   • d'éviter le cas d'une racine complexe avec un test si $x < 0$ qui va soulèver une erreur de valeur (ValueError)

   • de réaliser en fin de boucle, une comparaison avec l'expression $y = (1+\sqrt{1+4x})/2$ et d'afficher les résultats avec 16 chiffres décimaux (exemple :)

     n = 23
     r = 2.3027756377319943
     y = 2.3027756377319948
     

2. Examiner le comportement du code pour $x$ très petit ou nul ($x = 0$, $x \approx 0$)

3. Généraliser le code pour des racines nième ($n \ge 2$):

   $r_n = \sqrt[n]{x + \sqrt[n]{x + \sqrt[n]{x\ +\ ... }}} $

4. Discuter le cas ou $n<0$

In [26]:

#%% r = sqrt(x + sqrt(x + sqrt(x + ..))...)

x = 5

# test si x < 0 on soulève une erreur de valeur
if x < 0: raise ValueError("x doit être positif")

eps = 1e-8         # on fixe la précision

# on ajoute une sécurité pour éviter la boucle infinie
# en cas de non convergence
k , kmax = 0, 10000

s = x**0.5            # 1ere racine (terme 1)
r = (x + s)**0.5      # 2eme racine et plus (termes 2...n)

while (abs(r-s) > eps) and k < kmax:
    s = r
    r = (x + s)**0.5
    k += 1

print(f" n = 2")
print(f"{x = }")
print(f"{k = }")
print(f"{r = :.16f}")


## Validation : non demandée
y = (1 + (1+4*x)**0.5)/2
err = abs(1 - r/y)
print('\nvalidation')
print(f"{y = :.16f}")
print(f"{err = :.3g}")

 n = 2
x = 5
k = 11
r = 2.7912878468582938

validation
y = 2.7912878474779199
err = 2.22e-10

1.9. Exo 9 Hauteur d'eau - formule de Manning¶

Reprendre l'exemple 2.2 du cours ( Tirant d'eau en écoulement à surface libre - formule de Manning) avec la méthode de Newton-Raphson

On aura besoin de dériver : $$F(h) =\frac{1}{n} A(h) R(h)^{2/3} \sqrt{S} - Q.$$

soit : $$F'(h) =\frac{1}{n}\sqrt{S} \left(A'(h) R(h)^{2/3} + \frac{2}{3} A(h) R'(h)/ R(h)^{1/3} \right) $$

avec $A'(h) = b$ et $R'(h) = b^2/(b+2h)^2$

On peut alors soit coder ces trois expressions, soit les rassembler en une seule en remplaçant $A(h)$ et $R(h)$,et $A'(h)$ et $R'(h)$ dans $F(h)$ et $F'(h)$ :

$$ F(h)= \frac{\sqrt{S} b^{\frac{5}{3}} h^{\frac{5}{3}}}{n \left(b + 2 h\right)^{\frac{2}{3}}} - Q $$

$$ F'(h) = \frac{\sqrt{S} b^{\frac{5}{3}} h^{\frac{2}{3}} \left(5 b + 6 h\right)}{3 n \left(b + 2 h\right)^{\frac{5}{3}}} $$

In [27]:

# Paramètres du canal
# Méthode sécante à modifier en Newton-Raphson
b = 1.0       # largeur (m)
Q = 0.5       # débit (m3/s)
S = 0.0005    # pente
n = 0.015     # coefficient Manning
g = 9.81      # gravité

tol = 1e-10    # Tolérance
itmax = 20   # nombre max d'itérations

# initiaisation h = Hauteur critique
h = (Q**2 / (g * b**2))**(1/3)
A = b * h
R = (b * h) / (b + 2*h)
F = A * (1/n) * R**(2/3) * (S**0.5) - Q

print("Points initial : h critique ")
print(f"{h = :.6f}, {F = :.3e}")

it = 0           # initialisation du compteur de boucle

# Boucle while basée sur la tolérance et le garde-fou
while abs(F) > tol and it < itmax:

    # Calcul de F(h)
    A = b * h
    R = (b * h) / (b + 2 * h)
    F = A * (1/n) * R**(2/3) * (S**0.5) - Q

    # Calcul de dF
    dA = b
    dR = (b/(b+2*h))**2
    dF = (1/n) * (S**0.5) * (dA * R**(2/3) + (2/3) * A * dR / (R**(1/3)))

    h = h - F/dF

    it += 1
    print(f"It {it}: h = {h:.8f}, F(h) = {F:.3e}")

if it>=itmax:
    print("Nombre d'iterations maximal atteint : pas de convergence")
else:
    print(f"\nTirant d'eau par la méthode de Newton-Raphson : h = {h:.12f} m")

Points initial : h critique 
h = 0.294277, F = -3.574e-01
It 1: h = 0.81402217, F(h) = -3.574e-01
It 2: h = 0.74912958, F(h) = 5.552e-02
It 3: h = 0.74877523, F(h) = 2.998e-04
It 4: h = 0.74877522, F(h) = 9.967e-09
It 5: h = 0.74877522, F(h) = 1.110e-16

Tirant d'eau par la méthode de Newton-Raphson : h = 0.748775221110 m